The efficient computation of transition state resonances and reaction rates from a 

quantum normal form 
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A quantum version of a recent formulation of transition state theory in phase space is presented. 
The theory developed provides an algorithm to compute quantum reaction rates and the associated 
Gamov-Siegert resonances with very high accuracy. The algorithm is especially efficient for multi- 
degree-of-freedom systems where other approaches are no longer feasible. 
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Introduction. — The question of how, as Marcus for- 
mulates it, a system "skis the reaction slope" is one of 
the crucial questions in reaction dynamics. Experimental 
techniques like photodissociation of jet-cooled molecules, 
molecular beam experiments or transition state spec- 
troscopy give detailed information about the reaction 
process as has recently been demonstrated, e.g., for the 
'paradigm' reaction of Hydrogen atom-diatom collisions 
(see, e.g., the review paper 2]). A chemical reaction can 
often be viewed as the scattering problem across a sad- 
dle point of the interaction potential. The cumulative 
reaction probability is then given by 
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where t is the transmission subblock of the scattering op- 
erator S for energy E and the summation in the latter 
expression runs over all incoming reactant states with 
quantum numbers and outgoing product states with 
quantum numbers rip. The ab initio quantum mechan- 
ical computation of N{E) soon becomes very expensive 
if the number of atoms in the system increases beyond 
3 and one has to resort to suitable approximations. The 
main approach to compute N{E) classically is transition 
state theory which was invented by Eyring, Polanyi and 
Wigner in the 30's. The main idea is to define a di- 
viding surface that divides the energy surface into a re- 
actant and a product component and compute the rate 
from the directional phase space flux through this sur- 
face. In order not to overestimate the rate the dividing 
surface must not be recrossed by reactive trajectories. 
In the 70's Pechukas, PoUak and others Q showed that 
for two degrees of freedom such a dividing surface can 
be constructed from a periodic orbit that leads to the 
so called periodic orbit dividing surface. Recently it has 
been shown that a generalization to higher dimension can 
be achieved from a normally hyperbolic invariant man- 
ifold (NHIM) Q — a fundamentally new object that 
takes the place of the periodic orbit. The dynamics is 
controlled by the NHIM's stable and unstable manifolds 
which act as separatrices that divide the reactive trajec- 
tories from the nonreactive trajectories. The NHIM, its 



stable and unstable manifolds and a dividing surface with 
the desired properties can be directly constructed from 
an algorithm based on a Poincare-Birkhoff normal form 
procedure |5l|. 

Much effort has been devoted to developing a quan- 
tum version of transition state theory whose implementa- 
tion remains feasible for multi-dimensional systems (see 
the flux-flux autocorrelation function formalism by Miller 
and coworkers 0). In this Letter we present a quantum 
version of the normal form procedure that lead to the 
construction of the high-dimensional phase space struc- 
tures that govern the classical reaction dynamics and 
demonstrate that this quantum normal form approach to 
transition state theory provides an efficient procedure to 
compute quantum reaction rates and the corresponding 
Gamov-Siegert resonances . 

The quantum normal form. — The main idea of which 
the seed can already be found, e.g., in is to derive 
a local approximation of the Hamilton operator of the 
scattering problem that is valid near the saddle, and in 
order to facilitate further computations, takes a much 
simpler form than the original Hamiltonian. This can be 
achieved in a systematic way by a procedure based on 
the Wigner- Weyl calculus that has been used by others 
before to compute energy spectra associated with stable 
equilibria 9]. Here the manipulations of an operator A 
are expressed in terms of its symbol which is the function 
A(p, x) on the /-degrees-of-freedom phase space (p, x) G 
R2/ defined by 

(2) 

Defining the multiplication * of two symbols A and B 
according to 



(3) 



A*B^Aciip[i-[{dx,dp)-{dx,dp)] \B. 



where the arrows indicate whether the partial differenti- 
ation acts to the left (on A) or to the right (on B), gives 
the property that the quantization of the product of two 
symbols A and B is equal to the product of the quanti- 
zations of the individual symbols, i.e. AB = A* B . The 
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^-product leads to the definition of the Moyal bracket 

{A,B}m^^{A^B-B^A). (4) 

We define the order s of a monomial p°'x^K^ = 
Pi^ ■ ■ ■p'^^Xi^ ■ ■ ■ x'J^h" according to s = |a| + |/3| +2n = 
ai + • • • + a/ + /3i + • • • + /?/ + 2n, and denote the vector 
space of polynomials spanned by monomials of order s 
by W. For A e W and B e W' we define the Moyal 
adjoint by Mad^i B :— {A, B}m, then its iterates satisfy 
[MadA]"S G >v«('^-2)+'5'. 

We start from a Hamilton operator H whose symbol 
has expansion 

oc 

H = Eo + J2Hs (5) 

s=2 

where Eq is a constant energy and Hs € VV"*. Like in 
the case of classical transition state theory |5| we assume 
that the second order term is of the form 

i?2 = + + • • • (6) 
with the integrals 

I = Pix,, J, = l{pl+xl), k^2,...J. (7) 

This corresponds to a classical equilibrium point 
of saddle-center-. .. -center stability type ('saddle' for 
short), i.e. the matrix associated with the linearized 
system has one pair of real eigenvalues ±A associated 
with the saddle or 'reaction coordinate' and / — 1 pairs 
of imaginary eigenvalues iiw/t, k = 2, . . . , /, associated 
with the center or 'bath' degrees of freedom. Note that 
PiXi = (pi —ql)/2 where (pi,xi) and are related 

by a rotation of 45°. We restrict ourselves to the generic 
non-resonant case where the frequencies Wfc, fc = 2, . . . , /, 
are rationally independent 

In order to simplify the Hamiltonian H we will trans- 
form it by successive conjugations with unitary operators, 
H =: ^ H(3) ^ H(i) ^ , HiN)^ ^here 

fj(n) ^ gW„/r.j^(„-i)g-iw„/n ^ 

with Wn € W". Using the Moyal adjoint the symbol of 
the right hand side can be expanded as 

00 

^(n) ^ Y -M&dw I'^iff"-!) . (9) 

k=0 

If we expand furthermore each of the symbols i?'") in 
a generalized Taylor series as in ((SJ, i/^"^ = £'0 + 
J2s>2 with ifi"^ G then using ® the terms in 

these series can be related by 

^i"'- E fci[Madw-J'=i/l'!;;t,) , (10) 

k=0 



where [•] denotes the integer part. Notice that for s < n, 
Hs^'^ = Hi" , and for s = n we obtain 

Hi"^ ^H^r^^ +{Wn,H2} (11) 

(n) 

where we have used that H2 = H2 for all ?i > 3 and 
that the Moyal bracket reduces to the Poisson bracket 
{•, •} if one of the functions is quadratic. This is the 
homological equation which is familiar from the classical 
normal form algorithm, j5j] , and under the non-resonance 
conditions on the frequencies Wfe, given any 77^""^^ G W" 
there exists a unique W„ G W" such that i?^"^ can be 
written as a function of the actions J, J2, • • • , J/ alone 
(or, equivalently, satisfies {/, i/,^"'} — {J2,i?i"'} = 
••• = {J/,i?^)} = 0). 

Choosing the generators of the unitary transformations 
Wn recursively for n = 3, 4, ... as solutions of ((TT)) we 
obtain an operator H'^^^ whose symbol is of the form 

(^) ^Eo + Sf^2 + where the first part is 

a polynomial in the actions /, J2, • • • J/, i.e., is in normal 
form, and the remainder R^^^^"^ consists of terms of or- 
der -|- 1 and higher. In a final step we want to express 
the quantization of the normal form part as an operator 
function of the quantized actions / and Jfc. To this end 
we use a recursion relation for /" (and a similar one for 

the Jk), = i7" - /"-in2fi2/4, which can be derived 
from the product formula This allows us to express 
quantizations of powers of / as a polynomial in powers 
of /. In this way we find a polynomial ^^qnp such that 

U*^HUn = H^f^Fil -^2, • • • , Jf) + , (12) 

where Un = Ilia e"'"^"/'- ^qnf(^' ^2, - • , J/) is 
called the quantum normal form (QNF) of H of order 
N. The remainder term i?'-^'*'^-' has a symbol which is of 
order + 1 and is therefore very small near the saddle 
point. Hence the dynamics near the saddle point can be 
described with high accuracy by the QNF Hamiltonian. 
The advantage of the QNF Hamiltonian is that it is an 
operator function of the commuting operators / and Jk 
whose properties are well understood. Finally we point 
out that in the limit h we recover the classical normal 
form of order N, limri^o H^^p{I, J2, • • • , J/) = H^^\ 

We implemented the algorithm of the QNF computa- 
tion in the programming language CH — h. 

Resonances and reaction rates. — The eigenfunctions of 
the truncated QNF Hamiltonian are tensor products of 
harmonic oscillator wave functions for the center degrees 
of freedom Xk, k = 2, . . . , /, and eigenfunctions of the 
operator 

?=y(^i5.,+^) (13) 
associated with the saddle direction. The operator / has 
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FIG. 1: Transmission probability T{E) and Gamov-Siegert reso- 
nances in the complex energy plane for the Eckart potential. The 
parameters are a = 1, 10 A = B = 5 for the Eckart potential, and 
m = 1 and h= 0.1. 



eigenfunctions [lO| 

V'out±(xi) = e(±.Ti)|.Tirl/2+iVr. ^ (14) 

Q being the step function, which arc outgoing waves. In- 
coming waves can be defined from the Fourier transforms 



(2/i)e' 



xivi 



(15) 



Note that in order to classify the eigenfunctions as 'out- 
going' or 'incoming' it is convenient to rotate the coordi- 
nates {pi , xi ) back to the more standard notation (pi , qi ) 
for a potential barrier mentioned above. Expressing the 
functions V'ini in terms of 4'out± gives the entries of a 
'local S-matrix', 

'(pin+ = Sn 11 i'out ++ Sn 12 ^ont- , (16) 
Ipin- = Sn2li'out+ + Sn22^ont- ■ (17) 

Here n denotes the vector of nonnegative integers 
{n2,-..,nf) formed by the quantum numbers of the 
modes in the center directions. The local S-matrix is 
block diagonal due to the separability of the QNF. Mode 
mixing is a 'global' effect which occurs from connecting 
the local wave functions to the asymptotic reactants and 
products wave functions. However, the local S-matrix 
alone already contains the full information needed to 
compute reaction rates and resonances. 
Evaluating the integrals ifT^ gives 



e'(f--Ki"'^) /l / 



with / being implicitly defined by 

h{n2 + 1/2), ...,h{nf + 1/2)) = E . 
The transmission probability of mode n is 

r„(i?) = |5„i2(i?)p-(l + e-2-l)-i, 




(18) 
(19) 

(20) 
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FIG. 2: (a) Errors for the transmission probability of the Eckart 
potential computed from the QNF. (b) Errors for the QNF reso- 
nances (differences between the corresponding QNF and exact com- 
plex energies) of the Eckart potential as a function of the quantum 
number ni. The different colors correspond to different orders of 
the QNF. The parameters for the Eckart potential are the same as 
in Fig. [3 



which gives the cumulative reaction probability N(E) = 
J2n Tn{E). The S-matrix has poles at / = -\h{ni + 1/2) 
for nonnegative integers ni and these define the Gamov- 
Siegert resonances via (|19|) . 

Examples. — We at first illustrate the procedure for 
ID potential barriers, i.e. Hamiltonians of the form H — 
p^ I {2m) + V{x) , where V{x) has a maximum which we 
can assume to be at x = 0. The second order QNF is 
easily obtained and gives the well know result i/qNF, 2 — 
V{0) + A/ with A = (-F"(0)/m)i/2 which is equivalent 
to approximating the potential barrier by an inverted 
parabola. The first nontrivial correction to this result 
comes from the fourth order QNF i?QNF,4 given by 



V^(0) + A/ 
1 



1 



V"'{Q) + V^^{Q)\P 



IGm^A^ V BmA^ 
7 



T/"'(0) + y^^(0) ft- 



We apply the QNF to the Eckart potential 0| 



l^Eckart (x) = A 



B- 



\^Q(x+xo)/a (1 _|_ Q(2;-f2;o)/a'J 



(21) 



(22) 



with xo = a\n{B + A)/{B-A) and B > A>0. FigureHl 
shows the exact transmission probability which is known 
analytically, and the exact string of resonances together 
with the resonances from the second order QNF which 
have constant imaginary part. The bending of the string 
of exact resonances is a nonlinear effect that is very well 
described already by the fourth order QNF. The excellent 
accuracy of the resonances and the cumulative reaction 
probability computed from higher orders of the QNF is 
illustrated in Fig. [5] 

We next consider the two-degree-of-freedom example 
of a Hamiltonian with an Eckart potential in the xi- 
direction, a Morse potential 



VMorso(a;2) = De (c 



--2ax2 
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(23) 
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FIG. 3: (a) The top panel shows the cumulative reaction probabil- 
ity N(E) (oscillatory curve) and the classical flux divided by 27r?i 
(smooth curve) for the Eckart-Morse potential defined in the text 
with e = 0. The bottom panel shows the resonances in the complex 
energy plane marked by circles for the uncoupled case e = and 
by crosses for the strongly coupled case e = 0.3. For the coupled 
case the numerically exac t resonances are computed from the com- 
plex dilation method jj^l ■ The parameters for the Eckart potential 
are the same as in Fig. The parameters for the Morse potential 
are De = 1 a = 1. Again we choose m = 1 and h = 0.1. (b) 
Errors for the cumulative reaction probability in (a) for different 
orders of the QNF. (c) Average errors < |i?QNF — Eoxactl >ni = 
^,1.2=0 l-^QNp("-i,'^2) - i?cxact(ni, "2)1/6 for the QNF resonances 
in (a) for the coupled case e = 0.3 as a function of the quantum 
numbers ni. 



in the X2-direction, plus a kinetic coupling epiP2- In 
Fig. 131 the cumulative reaction probability and the res- 
onances computed from the QNF are compared with the 
exact results. In the uncoupled case e = 0, N{E) in- 
creases as a function of E at integer steps each time a 
new transition channel opens, i.e. when the transmission 
probability (E) of a Morse oscillator mode ri2 switches 
from to 1. For both the uncoupled and strongly coupled 
case the resonances form a distorted lattice parametrized 
by the mode quantum number n2 in horizontal direction 
and the quantum number ni in vertical direction. Similar 
to the ID example, each string of constant n2 is related 
to one step of N{E). 

Like in the ID case the agreement of the QNF results 
with the exact results is excellent and this remains the 
case even for the strongly coupled system. The QNF 
is an asymptotic expansion which in general does not 
converge. For the system shown in the example this can 
be seen from the fact that the 10th order QNF does not 
lead to an improvement of the results obtained from the 
8th order QNF. 

Conclusions. — The QNF computation of reaction 
probabilities and the corresponding Gamov-Siegert res- 
onances is highly promising since it opens the way to 
study high dimensional systems for which other tech- 
niques based on the ab initio solution of the quantum 



scattering problem like the complex dilation method 
or the utilization of an absorbing potential jl4| are no 
longer feasible. In fact, the numerical effort for com- 
puting the QNF is only slightly higher than the effort 
for computing the classical normal form (NF). The main 
difference is that the Poisson bracket in the classical NF 
needs to be replaced by the Moyal bracket. The storage 
of the QNF polynomials also requires only a disk space 
similar to the classical NF. Moreover, the QNF gives 
an explicit formula for the resonances from which they 
can be computed directly by inserting the corresponding 
quantum numbers. This leads to a direct assignment of 
the resonances. The QNF provides a quantum version 
of transition state theory that, in the semiclassical limit, 
is in accord with the classical phase structures that gov- 
ern the reaction dynamics. In fact, the classical phase 
space structures form the skeleton for the scattering and 
resonance wavefunctions, and exploiting this relationship 
which will give a deep insight into the now experimen- 
tally accessible quantum reaction dynamics is the subject 
of our future studies (see also 

This work was supported by the Office of Naval Re- 
search, EPSRC and the Royal Society. 



[1] 
[2] 

[3] 
[4] 

[5] 



[6] 
[7] 



[8] 

[9] 
[10] 

[11] 

[12] 
[13] 

[14] 



R. A. Marcus, Science 256, 1523 (1992). 

R. T. Skodje and X. Yang, Int. Rev. Phys. Chem. 23, 

253 (2004). 

P. Pechukas and F. J. McLafferty, J. Chem. Phys. 58, 
1622 (1973); P. Pechukas and E. Pollak, J. Chem. Phys. 
69, 1218 (1978). 

S. Wiggins, Physica D 44, 471 (1990); S. Wiggins, Nor- 
mally Hyperbolic Invariant Manifolds in Dynamical Sys- 
tems (Springer, Berlin, 1994); S. Wiggins, L. Wiesenfeld, 
C. Jaffe, and T. Uzer, Phys. Rev. Lett. 86, 5478 (2001). 
T. Uzer, C. Jaffe, J. Palacian, P. Yanguas, and S. Wig- 
gins, 957 (2001); H. Waalkens, A. Burbanks, and S. Wig- 
gins, J. Chem. Phys. 121, 6207 (2004). 
W. H. Miller, J. Phys. Chem. A 102, 793 (1998). 
R. S. Friedman and D. G. Truhlar, Chem. Phys. Lett. 
183, 539 (1991); T. Seideman and W. H. Miller, J. Chem. 
Phys. 95, 1768 (1991). 

R. Hernandez and W. H. Miller, Chem. Phys. Lett. 214, 
129 (1993); R. Hernandez, J. Chem. Phys. 101, 9534 
(1994). 

L. E. Fried and G. S. Ezra, J. Chem. Phys. 92, 3144 
(1988); P. Crehan, J. Phys. A 23, 5815 (1990). 
Y. Colin de Verdiere and B. Parisse, Comm. PDE 
19, 1535 (1994); Ann. Inst. Henri Poincare (Physique 
Theorique) 61, 347 (1994); Comm. Math. Phys. 205, 
459 (1999). 

S. Nonnenmacher and A. Voros, J. Phys. A:Math. Gen. 
30, 295 (1997). 

C. Eckart, Phys. Rev. 35, 1303 (1930). 
N. Moiseyev, Phys. Rep. 302, 211 (1998); W. P. Rein- 
hardt, Ann. Rev. Phys. Chem. 33, 223 (1982); B. Simon, 
Phys. Lett. A 71, 211 (1979). 

A. Neumaier and V. A. Mandelshtam, Phys. Rev. Lett. 



5 



86, 5031 (2001). (2005) 
[15] S. C. Creagh, Nonlinearity 17, 1261 (2004); ibid. 18, 2089 



